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Abstract: 



Relying on the recently proposed multicanonical algorithm, we present a numerical sim- 
ulation of the first order phase transition in the 2d 10-state Potts model on lattices up to 
sizes 100 x 100. It is demonstrated that the new algorithm lacks an exponentially fast in- 
crease of the tunneling time between metastable states as a function of the linear size L of 
the system. Instead, the tunneling time diverges approximately proportional to l? m . Thus 
the computational effort as counted per degree of freedom for generating an independent 
configuration in the unstable region of the model rises proportional to V 2 ' 3 , where V is the 
volume of the system. On our largest lattice we gain more than two orders of magnitude 
as compared to a standard heat bath algorithm. As a first physical application we report a 
high precision computation of the interfacial tension. 
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Critical slowing down is of crucial importance to computer simulations of phase tran- 
sitions. For second order phase transitions long autocorrelation times at criticality cause 
severe restrictions on the maximum lattice size for which one can obtain good statistics of 
thermodynamic quantities. For a number of spin systems this critical slowing down was 
overcome by the nonlocal cluster algorithm of Swendsen-Wang |IJ, for a recent review see 
||. However, for first order transition one encounters an even worse and different problem 
of critical slowing down. The interfacial free energy between disordered and ordered states 
has a finite value on the critical point for the infinite volume system. Configurations dom- 
inated by the presence of the interface will be exponentially suppressed by the Boltzmann 
factor in the canonical ensemble. On finite lattices this leads then to an exponentially fast 
suppression of the tunneling between metastable states of the system with increasing lattice 
size. To overcome this critical slowing down effect for first order transitions, we recently 
proposed a multicanonical Monte Carlo algorithm |§. The multicanonical MC algorithm is 
designed to enhance configurations, which are dominated by the presence of the interface 
and therefore exponentially suppressed. In this way it is possible to avoid the exponentially 
fast growing slowing down at the first order phase transition. In this paper we demonstrate 
this in the case of our example: the 2c? 10-state Potts model. 

The 2d 10-state Potts model |4j] is defined by the partition function 

z(p) = £ expoas), (i) 

configurations 

s = E <W ( 2 ) 

<i,j> 

q l ,q j = 0, ...,9. (3) 

Recently there has been renewed interest in this model || || [jj. It serves as an excellent 
laboratory for finite size scaling (FSS) studies of temperature driven strong first order phase 
transitions. We have chosen this model as the first testing ground for our new method by, 
essentially, two reasons: a) it has a strong first order transition and this is the situation for 
which our method promises the most dramatic improvements, and b) accurate data, due to 
recent large scale simulations |5|, [7j , exist in the literature. Beyond these reasons, an accurate 
determination of the interfacial free energy is certainly also of physical interest. 

To calculate the interfacial free energy F s between the disordered and the (ten) ordered 
states has remained the hardest problem. The reason is the pronounced double peak struc- 
ture of the sampled action density Pl{S) in the canonical ensemble near the critical point, 
as illustrated in Fig. 1 for lattices of size V = L 2 with L = 24 and L = 100. The figure relies 
on data obtained with the multicanonical MC algorithm to be discussed in this paper, and 
is arranged to correspond to the action density of the canonical ensemble at (3 = /3£. The 
pseudocritical point /3£ is defined such that both maxima are of equal height: 

P,l,max J-, I rYl,max\ j-> / rf2,raa:r\ r-)2,max / A \ 

L = ^L{b L ) = r L {b L ) = f L . (4J 

In addition we have imposed the normalization condition 1 = p^ max — p2,max p.^ ^ depicts 
the action densities for lattices with L = 16, 24, 34, 50, 70 and L = 100 on a logarithmic scale 
and we see that 4 orders of magnitude are involved: p™ m /pm ax ~ 5.1 X 10~ 5 for L = 100. 
With our conventions for Pl{S) the interfacial free energy F s = can now be defined || 
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as the L — > oo limit of the quantity 

F[ = —InPr. (5) 

For a numerical calculation of F[ it is now clear that any algorithm which samples config- 
urations with a probability ~ Pl{S) would slow down proportional to 1/P™ in . As for large 
lattices P|" m ~ exp(— F s L d ~ l ), it is expected that an appropriately defined tunneling time 
t\ will behave as 

T t L = A T L a e +FSLd '\ (6) 

The parameters A T and a can in principal be determined by a fit to the measured tunneling 
times. 

The multicanonical MC algorithm samples configurations with the weight 

Vr{S) ~ e { < + ^ s) for S k L <S < S k L +1 , (7) 

instead of sampling with the usual Boltzmann factor Pl(S) ~ exp(/?£S') corresponding to 
the canonical ensemble. Here we partitioned the total action interval < S < 2V into 
k = 0,...,N (N odd) intervals = (S'l, S k+l ]. The idea of the multicanonical MC algo- 
rithm is to choose intervals Ik and values of ft\ and a\ at the pseudocritical point /?£ in 
such a way, that the resulting multicanonical action density Vl{S) has a approximately flat 
behavior for values of the action in the interval [S l I f nax , s^ max ], that is to say: configura- 
tions dominated by the interface are no longer exponentially suppressed as they are in the 
metastable-unstable region of the canonical ensemble. Physically this can be achieved by 
choosing the /3-parameters fl\ such, that the system gets heated, when its in the ordered 
state of the metastable region, cooled when its in the disordered state, and neither of both if 
its in the unstable region. The parameters fil hereby take the form fil = ftl + Sftl, where the 
coupling constant difference 5 (3*1 changes sign as a function of S and is responsible for the 
altered dynamics of the model. The parameters a\ are adjusted in such a way that Vl{S) 
is a steady function of S. 

In we demonstrated, that when the double peak distribution Pl(S) can be approxi- 
mated by a double gaussian, the multicanonical action density Vl{S) can be made arbitrarily 
flat by driving a control parameter r > 1 towards 1. In this case we choose action values 
S k L with Si = 0, S? +1 = 2V } Si = Sl max , Sg = Sl max and in the interval [S]; max \ Sr n ) 
action values defined by the equation P L (Sl) = r 1 ^ k P L {S]f aax ) for k = 1, N/2. An analog 
procedure is adopted in the interval (5'£ im , S 2 £ max ]. Having defined the action values Sk and 
corresponding intervals Ik the setting 

f (3 C L for k = 0, N/2, N 
Pl= + Hr)/(S k L +1 - SI) for k = 1, N/2 - 1 (8) 
[ PI - ln(r)/(S£ +1 - Si) for k = N/2 + 1, N -1 

and the recursion 

= a k L+(P k L -P k L +1 )S k L + \ a° L = (9) 

defines the multicanonical ensemble. In accordance with detailed balance, standard Metropo- 
lis and heat bath updating algorithms have been generalized to the multicanonical situation 
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[§, |§. Finally we obtain the canonical action density distribution Pl(S) through a reweight- 
ing step similar to [[II], |IT| from the multicanonical distribution Vl{S): 

P L (S) = e^i-^-^VLiS) for S k L < S < S k L +1 . (10) 

As an example we show for our L = 70 system in Fig. 3 the multicanonical action 
density distribution V-m(S) together with the reweighted distribution P 7 q(S). In practice 
the appropriate choice of the parameters in eq. (7) is obtained by making from the given 
systems a FSS prediction of the density distribution Pl(S) for the next larger system. A 
second run may then be performed with optimized parameters. It is our experience that 
the guess works normally so well that the second run is barely an improvement as compared 
with the first. On the smallest systems standard MC simulation provides initial data. Our 
statistics for this investigation was 4 ■ 10 6 heat bath sweeps per run and lattice size. One 
sweep updates each spin of the lattice once. 

We define the tunneling time t 1 l as the average number of sweeps needed to get from a 
configuration with action S = S^' max to a configuration with S = S^ max and back. With our 
statistics of 4 ■ 10 6 sweeps per run the system tunnels then in total 8 ■ 10 6 /r| times. Table 1. 
collects the measured tunneling times. For our smaller systems we have also carried out 
standard heat bath MC runs at /?£ and the associated tunneling times are also reported in 
Table 1. For the larger systems standard MC runs would not tunnel often enough to allow for 
a reliable direct calculation of their tunneling times. This is of course due to the exponential 
slowing down of the standard MC simulation. In Fig. 4 we display on a log to log scale the 
divergence of the tunneling times r l L for the multicanonical MC algorithm (circles) and the 
heat bath algorithm (triangles). There is clearly a different behavior of the two algorithms 
involved. While for the multicanonical MC algorithm the increase of the tunneling time is 
consistent with a power law, the heat bath algorithm displays an exponentially fast growing 
tunneling time. Performing a x 2 — fit we obtain the following fits 

r| (multicanonical) = 0.73(3) ■ L 2 65{2) with x 2 /{n F - 1) = 0.96, (11) 
r* (heat bath) = 1.46 ■ L 2 ' 15 ■ e +a080xL with x7(™f - 1) = 1-34. (12) 

The quality of the fits as indicated by the x 2 values (rip — I denotes the number of degrees of 
freedom minus one) is reasonable. In case of the heat bath algorithm we could not reliably de- 
termine the errors from a 3 parameter fit. The ratio R = r^(heatbath) / ^(multicanonical) 
is a direct measure for the relative efficiency of the two algorithms. Using the fits we ex- 
trapolate its value to the L = 100 system and estimate a factor R « 500 for this case. The 
multicanonical algorithm approximately slows down like ~ y 2325 with respect to the number 
of updates per degree of freedom. This is only slightly worse than the optimal performance 
~ V 2 which was estimated in || based on a random walk picture. For the heat bath al- 
gorithm the inefficiency of the algorithm prohibits a very accurate estimate of F s from the 
behavior of the tunneling time according to eq. (6). The fitted value in eq. (11) is however 
close to the determination of the next paragraph (eq. (15)). 

Our multicanonical data allow the so far most precise determination of the interfacial 
free energy for the 2d 10-state Potts model. For this purpose we determine maxima and 
minima of the Pl{S) distributions by self-consistent straight line fits over suitable 5-ranges. 



3 



Together with the central values of their associated ranges, our F[ values are collected in 
Table 2. Performing the FSS fit according to M 

F[ = F s + ^ (13) 

we obtain consistent results for lattices of size L = 16, 24, 34, 50, 70 and L = 100, as displayed 
in Fig. 5 (we have x 2 / {np — 1) = 0.27). We estimate the infinite volume value of the interfacial 
free energy to be 

F s = 0.09822 ±0.00079. (14) 

This value may however depend weakly on the analytical form of the FSS fit || and even 
with our large lattices we may still face additional systematic errors of a similar order than 
the quoted statistical error. Future simulations on even larger lattice sizes might therefore 
be of interest. 

In summary, we have introduced a multicanonical ensemble for the numerical simulation 
of first order phase transitions, which eliminates an exponentially fast increase of the tun- 
neling time between the ordered and disordered states in the critical region of the system. 
This finding is achieved by replacing the usual equilibrium dynamics of the canonical en- 
semble, through a new equilibrium dynamics, where the ordered and disordered states of the 
system get heated and cooled in a well controlled way. Thus configurations dominated by 
the presence of the interface are enhanced during the simulation. 

The multicanonical MC algorithm gives a general framework for the numerical studies 
of first order phase transitions in statistical mechanics as well as for field theoretic models. 
From the numerical point of view the interesting question will be what improvement factors 
can be achieved as compared to standard algorithms for certain models on certain lattice 
sizes. We expect the answer to this question to be determined by the value of the quantity 
Q = F[ x where the strength of the first order phase transition is indicated by the 

magnitude of F[ and the d is the dimensionality of the system. In the case of the 2d 10-state 
potts model we find at values of Q ~ 10 approximately an improvement of two to three 
orders of magnitude, while at values of Q ~ 1 the improvement is marginal. 

An implementation of the multicanonical MC algorithm for non-Abelian gauge theories 
is straightforward and we think that future investigations of the QCD deconfining phase 
transition will benefit from this. Beyond first order phase transition, it may well be that 
multicanonical algorithms could be of use for other numerical calculations in statistical me- 
chanics, like estimates of the free energy or spin glass simulations. 
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Tables 



L 


t£ multicanonical 


t£ heat bath 


12 


542 (4) 


793 (7) 


12 


— 


776 (9) 


16 


1147 (10) 


1988 (23) 


24 


3354 (57) 


9634 (408) 


34 


8375 (245) 


43923 (3151) 


50 


23763 (1321) 


270565 (63222) 


50 


24932 (1064) 




70 


69492 (6383) 




70 


62218 (5560) 




100 


160334 (16252) 





Table 1: The tunneling times t 1 l as a function of the lattice size L for the multicanonical MC 
algorithm (second row) and the heat bath algorithm (third row). For some lattice sizes we 
display the results of several simulations, whose difference lies in slightly different coupling 
parameters. 



L 


Pi 


Ql,max 


Cjmin 


n'2,max 


n 


12 


1.40738 (09) 


116 


169 


243 


0.1071 (06) 


16 


1.41534 (12) 


216 


309 


429 


0.1086 (07) 


24 


1.42100 (08) 


523 


723 


978 


0.1058 (08) 


34 


1.42338 (09) 


1072 


1466 


1945 


0.1039 (13) 


50 


1.42481 (06) 


2358 


3162 


4192 


0.1027 (11) 


50 


1.42469 (06) 


2357 


3186 


4190 


0.1006 (10) 


70 


1.42536 (06) 


4661 


6257 


8190 


0.0983 (20) 


70 


1.42541 (05) 


4660 


6250 


8178 


0.1007 (12) 


100 


1.42576 (04) 


9602 


13060 


16686 


0.0986 (18) 


100 


1.42577 (04) 


9577 


13093 


16711 


0.0994 (15) 



Table 2: The pseudocritical couplings (5 C L , interfacial free energies and locations of the max- 
ima and minima of the action density distribution, as determined from the multicanonical 
distributions. 
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Figure captions 



Fig. 1 Canonical action density distributions Pl{s) for L = 24 and L = 100 lattices 
at the pseudocritical couplings /5|4 and (3^00. Here s denotes s — S/ (2V). The values of the 
maxima have been normalized to 1. 

Fig. 2 Action density distributions Pl{s) for lattices of size L = 16, 24, 34, 50, 70 and 
100 on a logarithmic scale. The values of the maxima have been normalized to I. 

Fig. 3 Multicanonical action density distribution V 70 (s) together with with its reweighted 
distribution P^s). The normalization of the distributions is chosen arbitrarily such that 
the figure looks nice. 

Fig. 4 Tunneling times for the multicanonical MC algorithm and the heat bath al- 
gorithm in a double log scale. The curves correspond to the fits in eq. (11) and eq. (12). 
The dashed part of the curve indicates the extrapolation to the L = 100 lattice for the heat 
bath algorithm. On the 100 lattice the systems still tunnels 50 times between the metastable 
states during 4 x 10 6 sweeps, when the multicanonical simulation is used. 

Fig. 5 FSS estimate of the interfacial free energy F s . Averages are used for those 
lattices for which we have two data sets. 
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